function [Z_samples,Pr,theta_samples,log_importance] = ...
        run_modQTL_association_pip(EXPR,GENO,ANNOT,CIS,NSAMPLE,THETA_MIN,THETA_MAX,NEPOCH)
    %
    % Estimate posterior inclusion probability (PIP) on eacn SNP
    % using importance sampling of hyper-parameters from
    %  theta ~ [THETA_MIN, THETA_MAX] uniformly.
    %
    % [Z,Pr,Theta,ELBO] = PIP(EXPR, GENO, ANNOT, CIS, NSAMPLE, THETA_MIN, THETA_MAX, NEPOCH)
    %
    % EXPR      = gene x individual expression
    % GENO      = snp x individual genotype
    % ANNOT     = snp x K prior SNP annotation
    % CIS       = gene x snp cis-relationship
    % NSAMPLE   = importance sampling time
    % THETA_MIN = min of theta value
    % THETA_MAX = max of theta value
    % NEPOCH    = optimization steps (with constraint)
    %
    % Z         = snp x samples conditioned PIP
    % Pr        = sample x 1 normalized probability
    % THETA     = sample x (K+1) prior log-odds parameters
    % ELBO      = sample x 1 (log-importance)
    %
    % code: Yongjin Park, ypp@csail.mit.edu
    %
    % Note: Importance sampling idea was borrowed from Carbonetto &
    % Stephens, Bayesian Analysis (2012)
    %

    % ****************************************************************
    [Ngene,Nind] = size(EXPR);
    [Nsnp,~] = size(GENO);

    ANNOT = [ones(Nsnp,1),ANNOT];
    K = size(ANNOT,2);

    TOL = 1e-2;
    SIGSQ = var(EXPR(:));
    width = abs(THETA_MAX - THETA_MIN);

    % ================================================================
    % common starting point

    theta_init = THETA_MAX*ones(K,1);
    prior_init = ANNOT * theta_init;

    fprintf(1,'\n\nFind common starting point\n');

    CIS = logical(CIS);

    [gene_idx, snp_idx] = find(CIS);
    [W_init,~,Z_init,invtausq_w_init,llik,elbo_init] = ...
        modQTL_association_data(EXPR,GENO,CIS,prior_init,SIGSQ,10,NEPOCH);


    % W_sum = 0*W_init;
    Z_samples = zeros(Nsnp,NSAMPLE);
    log_importance = zeros(NSAMPLE,1);
    theta_samples = zeros(K,NSAMPLE);

    % explore regular grids
    theta = repmat(THETA_MIN,K,1);

    for ii = 1:(NSAMPLE^K),

        theta_rand = THETA_MIN + width*rand(K,1);

        prior_odds = ANNOT * theta_rand;

        [~,~,Z,~,llik,elbo] = modQTL_association_data(EXPR,GENO,CIS,prior_init,SIGSQ,0,NEPOCH,...
                                                      W_init, Z_init, invtausq_w_init);

        theta_samples(:,ii) = theta_rand;
        log_importance(ii) = elbo;
        Z_samples(:,ii) = Z;

        fprintf(1,'Importance Sampling %03d, LLIK = %.4e, ELBO = %.4e, Theta = [%.2e, %.2e], Pr{Z} = %.2f\n',...
                ii,llik(end),elbo,min(theta_rand),max(theta_rand),mean(Z));

    end

    Pr = exp(log_importance - max(log_importance));
    Pr = Pr / sum(Pr);

end
